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REGENERATIVE  SIMULATION  OF  HARRIS  RECURRENT  MARKOV  CHAINS 


Peter  W.  Glynn 


I .  Introduction 

In  recent  years,  a  substantial  amount  of  research  effort  has  been 
devoted  to  development  of  efficient  techniques  for  statistical 
analysis  of  simulation  output.  When  the  stochastic  system  being  simu¬ 
lated  enjoys  regenerative  structure,  the  regenerative  property  can  be 
exploited  to  produce  a  methodology  which  possesses  many  attractive 
properties.  Such  an  approach  has  proven  useful  in  the  steady-state 
simulation  problem  (CRANE  and  IGLEHART  (1974),  FISHMAN  (1973)), 
quantile  estimation  (IGLEHART  (1976)),  and  selection  of  the  "best" 
system  (IGLEHART  (1977)).  The  applicability  of  these  procedures  has 
been  limited,  however,  by  the  need  to  construct  regeneration  times  for 
the  process  under  study. 

In  this  chapter,  we  will  extend  the  regenerative  method  of 
simulation  output  analysis  to  the  class  of  Harris  recurrent  Markov 
chains.  These  chains  are  a  natural  setting  for  an  extension  of  this 
kind,  since  it  can  be  shown  that  every  generalized  semi-Markov  process 
(GSMP)  corresponds  to  a  certain  Markov  chain  taking  values  In  a 
complete,  separable  metric  space.  As  WHITT  (1980)  has  pointed  out, 
GSMPs  are  a  mathematical  formulization  of  the  general  discrete-event 
simulation.  It  is  shown  in  [12]  that  if  the  Markov  chain  correspond¬ 
ing  to  the  GSMP  is  well-behaved  from  a  simulation  standpoint,  then  the 
chain  must  in  fact  be  Harris  recurrent.  Thus,  the  class  of  Harris 
chains  contains  the  class  of  well-behaved  discrete-event  simulations. 
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In  Section  2,  we  will  briefly  describe  an  embedding  for  Harris 
chains  that  provides  a  weakly  regenerative  environment  (a  weakly 
regenerative  process  generalizes  the  concept  of  regenerative  process 
in  the  sense  that  the  sequence  of  tours  may  be  m-dependent) .  Section 
3  develops  an  algorithm  ,  very  similar  to  that  of  the  regenerative 
method,  for  statistical  analysis  of  the  steady-state  simulation 
problem.  Section  4  concerns  a  refinement  of  the  algorithm  for  the 
case  where  the  embedding  is,  in  fact,  regenerative. 

For  a  given  Harris  chain,  there  are,  in  general,  many  ways  of 
embedding  the  process  in  a  weakly  regenerative  environment.  In 
Section  5,  we  investigate  the  existence  of  “optimal"  embeddings. 
Section  6  concerns  application  of  the  ideas  of  Sections  2  through  5'  to 
the  "passage  time"  problem  for  Harris  chains.  Recently,  there  has 
been  substantial  interest  in  passage  time  problems  for  networks  of 
queues;  see,  for  example,  IGLEHART  and  SHEDLER  (1980).  We  consider 
this  concept  in  the  Markov  chain  setting,  and  prove  that  one  can 
reduce  the  passage  time  problem  for  Harris  chains  to  the  situation 
covered  in  Secions  2  through  5.  Finally,  in  Secion  7,  we  consider  a 
simple  vector-valued  storage  process  operating  in  discrete  time,  and 
prove  Harris  recurrence.  This  storage  process  provides  examples  of 
chains  which  can  be  made  regenerative  only  via  an  embedding  of  some 
kind  (such  as  that  developed  in  Sections  2  and  3).  Thus,  the 
techniques  considered  here  do  legitimately  extend  the  classical 
regenerative  method.  ■  •— y -  T 

K 
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2.  Harris  Recurrent  Markov  Chains 

Let  E  be  a  complete,  separable  metric  space  with  E  its 
associated  Borel  sets.  Given  a  probability  transition  kernel  P  and 
a  probability  p  on  (E ,JE ) ,  one  can  construct  a  probability  measure 
P^  on  the  product  space  Q=*ExEx**»  ,  equipped  with  product 
o-field  F»ExEx««.  ,  such  that 


(2.1) 


Vx0  c  V  X1  c  A1*  •••’  Xn  c  V 


-  /  p(dxQ)  /  P(xQ,  dXj)  ...  /  P(xn_r  dxn) 
A0  A1  An 


for  Af  c  _E,  where  X^Cdi)  **  is  the  projection  onto  the  ith 
coordinate  of  Q.  The  process  {Xn}  is  called  the  Markov  chain 
associated  with  kernel  P  and  initial  probability  p.  We  refer  the 
reader  to  OREY  (1971)  for  notation  and  additional  properties. 

The  kernel  P  is  said  to  be  a  Harris  kernel  if  there  exists  a 
probability  $  on  (E,j£),  a  positive  integer  m,  a  non-negative 
^-measurable  function  X(x) ,  and  a  set  A  e  E  such  that: 


(2.2) 


The  set 
denoted 


i)  Pm(x,.)  \(x)  $(•),  for  all  x  c  E, 

ii)  inf{X(x):  x  e  A}  -  A  >  0, 

iii)  P  {S,  <  *}  -  1  for  all  probabilities  p  on  (E,E) , 

P  A  — 

where  S  -  inf{n  >0:  X  c  A} . 

A  “*•  II 

iv)  P  (X  c  A  infinitely  often}  «  i  for  all  x  c  A. 
x  run 

of  all  pairs  ($,X)  satisfying  (2.2)  for  some  A  will  be 
by  Muj.  The  measure  $  will  be  referred  to  as  a  splitting 
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measure.  Finally,  any  chain  for  which  its  kernel  is  Harris  will  be 
called  Harris  recurrent. 

It  can  be  shown  that  the  Harris  condition  (2.2)  is  equivalent  to 
requiring  existence  of  a  a-finite  measure  v  such  that  v(A)  >  0 
implies 

(2.3)  P^{X  C  ^  infinitely  often}  *  1 

for  all  x  c  E  (see  Section  4,  (12]).  All  Harris  chains  Possess  a 
unique  o-finite  invariant  measure  x.  This  measure  has  the  property 
that  x(A)  >0  if  and  only  if  (2.3)  holds  for  every  x  c  E  (see 
Theorem  2.7,  REVUZ  (1975)).  We  will  use  these  equivalent  formulations 
in  our  analysis  of  the  “passage  time"  problem  In  Section  6. 

Harris  chains  are  of  great  importance  in  simulation.  As  WHITT 
(1980)  has  indicated,  a  reasonable  mathematical  formalization  of  the 
general  discrete-event  simulaiton  is  the  class  of  generalized  semi- 
Markov  processes.  As  is  well-known,  such  processes  can  be  regarded  as 
Markov  chains  taking  values  In  a  complete,  separable  metric  space  (see 
(25]  for  a  discussion).  On  the  other  hand,  any  chain  for  which  the 
associated  steady-state  simulation  problem  is  well-posed  must  neces¬ 
sarily  be  Harris  recurrent  (see  [12],  Proposition  3.8  for  definition 
and  results).  Hence,  the  study  of  well-behaved  discrete-event  simula¬ 
tions  leads  naturally  to  the  study  of  Harris  chains. 

Any  Harris  chain  can  be  embedded  in  a  weakly  regenerative 
environment  (see  Section  4  of  (12]  for  definitions  and  details).  The 
construction  of  such  an  environment  proceeds  as  follows. 
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Under  our  topological  assumptions  on  (E,E),  there  exist  regular 
conditional  distributions  Px(y,dz)  such  that 

'  VB>  “  J  [J  Px<y,dz)]  Pm(x,dy) 

for  all  B  e  E®,  the  m-fold  product  o-field  of  ^(Px( *)  Is  the 

probability  on  (Q,£)  associated  with  initial  distribution 

concentrated  at  x).  Assume  that  is  non-empty,  and  let 

(4>, A)  c  Mffl»  Put  *  Em  x  {0,1}  and  let  be  the  corresponding 

Borel  cr-field.  For  x  c  E  and  B  c  E™,  define  the  kernel  P  by 
the  formlula 

(2.4)  P(x,  B  x  {0}) 

•  (1  -  X(x))  /  (/  lg(v,u)  Px(u,dv)]  Q(x,du) 

P(x,  B  x  {1}) 

-  X(x) )  /  [/  IB(v,u)  P^(u  ,dv)  ]  ij»(du) 
where  Q(x,  du)  is  defined  by  the  relation 

Pm(x,-)'«  X(x)  $(•)  +  (1  -  X(x) )  Q(x, •)  . 

This  kernel  defines,  for  each  x  c  E,  a  probability  i?x  on 

E  xE  x  • • •  as  follows: 

m  m 
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(2.5)  Px{Y1  o  B L ,  ....  Yn  £  Bn,  61 


»  •  •  •  » 


-  i, 


6  -  i  } 
n  a 


/  P(x,  dy.  X  {i  })  •< 
B1 


/  ^yn-l(m)’  dyn  *  {in}) 

n 


Here,  (Yi,61)  is  the  ith  coordinate  projection  on  Em  *  Em  *  ***  * 
where  Y  is  the  m-vector  CY±Cl) »  ....  Y^m)).  Finally,  let 
E*  -  u*  o  E^  x  A»  A  as  in  (2.2),  and  let  P^  be  the 

probability  on  E*  x  Em  x  En  x  ...  given  by 


(2.6)  P  {(YQ(0) . Y0(k))  c  Bq;  N  -  k;  (Yj  ,  b^ ,  ....  Yn,  6Q)  £ 

"  /  PU^X0’  •••’  *k^  c  dy0*  SA  “ 


*  Py0(k){(Yl,6l,"',Yn'6n'>  C  V 


wt  »re  Yq  ■  (YQ(0),  .  . 
and  SA  ■  inf{k  0: 

{ Xn}  by  the  relation 


Yq(N>)  is  the  coordinate  projection  on  E*, 
<r  A) .  We  now  define  the  random  variables 


-  Y0(k);  k  <  N 

<2'7)  Wi  ‘  Vi«>; 

The  process  {X^}  records  the  consecutive  components  of  the  Y^'s. 

It  is  not  difficult  to  show  that  {Xn}  has  marginal  distribution 

P  on  (1,F)  and  hence  (2.7)  acts  as  an  embedding  of  (X  }  within 
\i  ~  n 

the  process  YQ ,  (Yj.6^,  ...  (see  Proposition  4.10,  [12)). 
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The  process  Yg ,  (Yj,6^),  ...  is  weakly  regenerative  of  order  1 
with  respect  to  the  random  times  {T^ :  j  _>  2},  where  Tq  *  0  and 

T,  -  inffk  >  Vl,  \  -  1)  ■ 

Basically,  this  means  that  the  random  "tours" 

((Tk+rTk’  (yt  *  5t  >»  ••••  <YT  -i*  &r  k>  2) 

rk  k  Tk+1  1  ~k+l  1 

are  identically  distributed  and  1-dependent  (i.e.,  tours  k  and  j 
are  independent  for  |k-j|  >  1);  see  Section  4  of  [12J  for  a  more 
precise  description  of  the  result.  One  can  therefore  regard  the  prob¬ 
ability  space  constructed  above  as  a  weak  regenerative  embedding  of 
the  chain  {Xn}. 

Before  turning  to  a  simulation  implementation  of  the  above 
embedding,  let  us  recall  that  we  can  decompose  Pm(x,*)  into  its 
^-singular  and  ^-absolutely  continuous  parts,  yielding 

Pm(x,B)  =  /  h(x,y)  $(dy)  +  p“(x,B) 

B  S 

where  P™(x,»)  is  the  ijj-singular  part  of  the  decomposition.  Of 
course,  by  (2.2)(i),  h(x,y)  >  l(x)  for  $  a.e.  y.  We  can  now  state 
the  simulation  implementation  of  the  above  weak  regenerative  embed¬ 
ding. 

A 

1.  Set  Tg+0,  k<-0,  j  f  1. 

2.  Generate  Xg  with  distribution  p.. 

3.  If  X^  e  A,  go  to  5. 
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Set  k  <-  k+1 


4.  Generate  X^+j  from  the  distribution  P(X^,»). 

Go  to  3. 

5.  Set  Z  Xk,  A  *■  k. 

6.  Generate  from  the  distribution  P(X^,*).  Set  k  «-  k+1. 

7.  If  k  <  A+m,  go  to  6. 

8.  If  KZ.Xjj)  -  0,  go  to  5. 

9.  Generate  a  random  variable  U  uniform  on  [0,1]. 

10.  If  \(Z)/h(Z,Xk)  £  U,  go  to  5. 

11.  Set  Tj  «■  1+1,  j  *  j+1. 


The  random  times  {T^}  produces  by  this  algorithm  are  weak 
regeneration  times  for  the  process  {Xn}.  It  should  be  pointed  out 
that  each  random  variable  in  the  above  algorithm  must  be  generated 
independently  of  all  other  random  variables  (modulo  the  fact  that  the 
distributions  involved  are  correlated). 

The  above  algorithm  incorporates  one  major  modification  to  the 
embedding  given  by  (2.3)  through  (2.7),  namely  the  use  of  “acceptance- 
rejection"  to  calculate  those  times  at  which  {Xn}  is  distributed 
according  to  $(•)•  Let  U  be  a  uniform  variate  on  [0,1],  and 
observe  that 


P  {(X.  ,  ... ,  X  )  e  B;  h(x,  X  )  >  0;  U  <  \(x)/h(x,  X  )} 

x  i  m  ID  —  ID 

-  ex(ia  px(u  <  K(x)/h(x,xm)|x1,  ....  Xm>> 

-  Ex(Ia  \(x)/h(x,Xffl)} 
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where 


a  -  {(X. ,  . .. ,  X  )  e  B;  h(x,  X  )  >  0} 

1  m  m 

The  above  expression  can  be  simplified  as 

VVA|V  X<x)/h(x.  xm>> 

-  E  { /  I„( v ,X  )  P  (X  ,dv)  \(x)/h(x,  X  ) ;  h(x,X)  >  0} 
x  omxm  m  m 

*  X(x)  /  (/  IgCv.u)  Px(u,dv>]  i{i(du)  , 

which  is,  we  observe,  half  of  relation  (2.3).  The  other  half  of  (2.3) 
is  equally  easy  to  verify,  and  thus  we  see  that  generation  of  the 
additional  variable  U  allows  us  to  avoid  explicit  generation  of 
random  vectors  with  distribution  Px(u,  dv)*  The  term  "acceptance- 
rejection"  derives  from  the  interpretation  of  U  as  being  a  variable 
with  indicates  when  to  "accept"  a  variate  as  having  come  from  the 
distribution  <t>(*). 

3.  Confidence  Interval  Estimation  for  the  Steady-State  Simulation 
Problem 

Every  Harris  recurrent  chain  (Xn)  with  invariant  probability 
u  has  the  property  that 

n 

(3.1)  l  f(X.  )/n  *  wf  =  /  f(y)  n(dy) 

k-0 
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Px  a.s.  for  all  x  c  E,  provided  x|f|  <  •  (see  Proposition  3.8  of 
[12]).  Given  that  (Xn)  Is  the  output  of  a  simulation,  a  simulator 
is  often  interested  in  obtaining  point  estimates  and  confidence 
intervals  for  the  parameter  xf;  this  is  known  as  the  steady-state 
simulation  problem.  We  will  now  show  that  the  weak  regenerative 
embedding  of  Section  2  can  be  used  to  advantage  in  obtaining  a 
solution  to  the  above  problem. 

The  strong  law  (3.1)  immediately  provides  a  strongly  consistent 
point  estimator  for  xf.  Confidence  interval  estimation  depends  on  a 
central  limit  theorem  (CLT)  for  the  summands 
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where 


— >  denotes  weak  convergence,  N(0,  o^)  ia  a  normal 
distribution  with  Bean  0  and  variance  o^,  and 

a]  -  (E  (Y -(f))2  +  2E  (Y  (f)  Y -(f)))  . 

1  p  2  p  2  3 

Proof.  The  sequence  (Y^(f);  1^2}  is  a  sequence  of  identically 

distributed  1-dependent  random  variables,  so  we  may  apply  Theorem  20.1 

of  BILLINGSLEY  (1968)  to  obtain  (3.3)  for  cr^  >  0.  On  the  other  hand, 
2 

if  Oj  »  0,  then  it  is  easily  computed  that  the  variance  of 

T  -1 

—1/9  ^  A 

(3.4)  n  i/Z(  l  fl\)) 

k-T2 

converges  to  0,  and  hence  (3.4)  converges  to  zero  in  probability.  # 

This  suggests  the  following  simulation  algorithm. 

1.  Simulate  {X^}  t0  time  Tn+2' 

2.  Compute  Y2(f),  X2>  •••,  Yn+1(f),  x^. 

3.  Compute 
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„  o+l  n+1 

-  l  Vf>/ 1  \ 

k-2  k-2 


,  n+1  .  tri-1  n  . 

Sn  "  l  \(f)/n  ~  2rn  I  \  \(f)/n  +  rQ  l  \/n 
n  k-2  n  k-2  k-2 


:n  *  J-  Yk(f)  Yk+l^f)/n  "  rn  J  \  Yk+l(f)/t 
k-2  k-2 


"  "n  I,  Yk(f)  Wn  +  "n  l  \Wn 

k»2  k*2 


v  »  s  +  2c 
n  n  n 


4,  Let  $(x)  ■  P{(N(0,1)  x}  ,  put  25  «  1-6/2),  and  set 


L 

n 


z6  Vn/2/,xn  n 


1/2 


R 

n 


r 

n 


1/2,-  1/2 
z,  v  /t  n 
6  n  n 


where 


n+1 


We  claim  that  the  random  interval  [Ln,  1^)  is  an  approxi- 
mate  100(l-£)£  confidence  interval  for  r(f),  provided  >  0.  This 
follows  from  Proposition  3.2  and  the  "converging  together"  lemma 
(Theorem  4.1,  [2]),  together  with  the  observation  that  Ex 2  <  ® 
when  {Xn}  possesses  an  invariant  probability  (see  Proposition 
4.15,  [12]). 
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The  above  simulation  algorithm  requires  that  the  time  horizon  of 
the  simulation  be  random,  namely  T  j  time  units.  The  following 
algorithm  deals  with  fixed  length  simulations. 

1.  Simulate  Xq ,  Xj ,  ...,  Xn. 

2.  Let  l( n)  *  max{k:  T,  <  n} ,  and  compute  t v„, 

k  —  r  i(n)  i(n) 

3.  Set 


T  *  1/2  ,-1/2  1/2 

n  “  ri(n)  Z5  Vi(n)  'i(n)  n 


1/2  ,-1/2  1 

ri(n)  +  Z5  Vl(n)/TJl(n)  n 


(if  i(n)  <  2,  report  L^  =  =  T”  f(X^)/n).  Again,  the 

interval  [Ln,Rn]  is  an  approximate  100(1-6)%  confidence 
interval  for  itf;  this  follows  from  the  asymptotic  validity  of 
the  first  interval,  application  of  the  random  time  change  theorem 
(Theorem  17.1,  [2]),  and  another  application  of  the  "converging 
together"  lemma. 


4.  Regenerative  Embeddings  for  a  Special  Class  of  Harris  Chains 
Let  {Xn}  be  a  Harris  chain  such  that  is  non-empty. 

Then,  several  simplifications  of  Sections  2  and  3  are  possible. 

Let  P  be  the  kernel,  corresponding  to  (<J>,\)  e  ,  defined  in 
Section  2,  and  let  (Xq,6q),  (Xj,6j),  ...  be  the  coordinate 
projections  on  x  x  ....  Then,  take  P^  as  the  probability  on 
Ei  x  Ei  x  •••  defined  by 
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VX0  C  B0’  **•’  Xn  C  Bn*  60  ”  i0 . 6n  "  V 


/  p(dxQ  x  {iQ})  /  P(xQ,  dxx  x  ...  /  P(xn_lt  dxn  x  {in} 


where  Bq,  Bj,  c  _E  and 


h(Bq  x  {0})  -  ii(Bq )  . 


It  is  easily  verified  that  {Xn}  has  marginal  distribution  P^. 
Furthermore,  if  we  put  Tq  m  0, 


ra  =  inf{k  >  Vl:  \  =  1}  , 


the  process  {Xn}  is  regenerative  with  respect  to  the  random  times 
{Tj}  (i.e.,  weakly  regenerative  of  order  0).  The  proof  of  this 
result  is  similar  to  that  of  Proposition  4.11  of  [12]. 

The  simulation  algorithm  of  Section  2  may  be  appropriately 
simplified: 

1.  Set  Tq  -*-0,  k  ■«-  0,  j  *■  1. 

2.  Generate  Xg  with  distribution  p. 

3.  Sec  Z  *■  Xjj. 

4.  Generate  from  the  distribution  P(X^,«).  Set  k  «-k+l. 

5.  If  h(Z,  Xfc)  -  0,  go  to  3. 

6.  Generate  a  random  variable  U  uniform  on  [0,1]. 

7.  If  \(Z)/h(Z,Xk)  <  U,  go  to  3. 

8.  Set  Tj  *■  k,  j  *■  j+1. 


Assume  for  the  remainder  of  this  section  that  the  invariant 


measure  n  is  a  probability.  To  obtain  a  confidence  interval  for 
itf,  we  consider  the  sequences 


V«  '  I  f(V 

k‘Ti 

•  IjCl)  . 

For  a  simulation  on  the  random  time  interval  {0,  1,  Tn+i>,  we 

proceed  as  follows: 

1.  Compute  Yj(f ) ,  x ^ ,  ...»  Y^Cf), 

2.  Compute 


:n  -  l  Vf)/ 1  Tk 

n  k-1  K  k-1  K 


,  n+1  .  n+1  n  _ 

8n  *  X  Vf)/n  -  2rn  X  \  Yk(f)/n  +  rn  X  Tk/n 

k-2  k-2  k-2 


Set 


L  *  r,/  >  -  z.s  /t  n 
n  i(n)  6  n  n 


1/2 


R  -  r„,  .  +  z.s  /r  n 
n  !(n)  6  n  n 


1/2 


where 
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Then,  [Ln,Rn]  is  an  approximate  100(1-6)%  confidence  interval  for 
ref,  provided  rt  |  f  J  <  ®  and  0  <  Y2(|f|)  <  “. 

For  a  fixed  length  simulation,  let  l( n)  -  max{k:  n},  and 

use  the  interval  [Ln,Rn]  given  by 


Ln  ri(n) 


Z6SJl(n)/^Jl(n) 


-1/2 


n1/2 


R 

n 


,-1/2  1/2 

ri(n)  +  Z6SJL(n)/xJL(n)  n 


The  proofs  of  thse  CLT's  follow  from  the  fact  that  the  constant  cr^ 
of  Proposition  3.2  reduces  to 

2  *  *  2 

°1-EP  Vf) 

in  the  regenerative  case. 

-  -  2 

Given  the  important  role  that  Y^Cf)  plays  in  the  central 
limit  problem,  it  is  convenient  to  reduce  the  question  of  finiteness 
of  this  parameter  to  one  involving  parameters  defined  on  the  original 
{ Xn)  process  alone. 


(4.1)  PROPOSITION.  Assume  that  is  non-empty  and  that  ($,\) 

and  K  are  aa  in  (2.2),  with  a  »  l.  Than,  CY2C  jf|)?  <  *. 

provided  that: 

1)  «|f|2r  <  ..  r  >  X 
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ii)  {Rx  log  Rj}28  <  -,  a  >  1 

where  Rj  -  lnf{n  >^1:  Xq  e  A> ,  *A(*)  “  x(*  n  A)/x(A),  and 

1/r  +  l/s  >  1. 


Proof.  By  condition  (2.2)(i)  and  (ii), 


<  /  P(x,.)  u.(dx) 


and  hence 


V1 


Tl-1 


EhU  lf(xk>|>  -  t  Jo  lf(xk>|}‘ 


Trl 


J  Ex{  J  |f(Xk)| y  ♦(dx) 


V1 


<  //  Ey{  I  |f(xk)|}  *A(dx)  P(x,dy)/X 


Tr1 


E,  lEx  I  I  |fUk)|}'  A 

*A  1  k-0  1  K  1 


Tj-1  T2-l 

■  \(  j,  i«v|)2a<\(  j0  i«m>2/x 


where  *  lnf{k  >_  2:  ^  “  1}.  Let  ■  0,  R^  ■  inf  (k  >  R^_ 

c  A}.  Clearly,  the  last  term  above  can  be  bounded  by 
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®  n-1  A  n+1 


C4-2)  I  l  E*  {(  l  |f(Xk)|V ;  &R  +1  -  1.  6  -  6. 

n-2  i-1  11 A  k-0  1  *  1  n  Rj 


•  n-1  * 


n-1 


<  l  l  \  U°+D  I  V  ;  6  -  6  : 

n-2  i-1  * A  1-0  1  j  J 


where  we  have  used  the  inequality 


j  <  n}/X 

<  n}/X 


n 

<  I 

i-1 


XiV 


n 

<n  I 

i-1 


x 


2 

i 


in  the  second  step  and  Vj^  denotes  the  quantity 


Vr1 


I  |f(Xk) 


k-R 


For  the  individual  term  in  the  above  iterated  sum,  we  have,  by 
Holder's  inequality,  a  bound  of 


(4.3) 


l  K  <v?p> U? 


i-o  V  1 


*,  <V»  ■  61  rJ< 

A  J 


where  1/p  +  1/q  »  1  and  p  >  1.  A  simple  "geometric  trial"  argument 
(see  Proposition  4.11,  [12]  for  an  example)  shows  that  we  can  bound 
the  6„  .  terms  by  (1-X)^n  on  the  other  hand,  if  we  use  the 

K  J  + 1 
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fact  that  is  an  invariant  probability  for  the  chain  {X^} 

(see  (20),  p.  32),  then  we  obtain 


K  {V5P>  *  K  {Ex  ^vnP)>>  “  <vnP>  • 

UA  1  nA  h  .  °  *A  ° 

n-1 

Combining  these  bounds  with  (4.3),  and  substituting  in  (4.2),  gives 
the  bound 


2  T  n3  E  {V2p}1/p  (l-K)(n_1)/q/\ 

Tt.  0 

n*2  A 

a  2 

so  that  "C Y£(  |  f  | ) }  <  ®,  provided  EuA^V0P^  ^  ”•  Duplicating  the 

reasoning  of  COGBURN  (1970),  p.  506,  shows  that  conditions  4.1(i), 

(ii)  suffice  to  guarantee  that  E  V2p  <  ®.  J 

*A  0 

Proposition  4.1  shows  that  if  f  is  a  bounded  function,  then  the 

24.5 

central  limit  theorem  holds  if  E  R7  <  »  for  some  6  >  0. 

*  A  1 

The  confidence  intervals  introduced  thus  far  in  this  section  make 
explicit  use  of  the  regenerative  structure  of  the  process  {Xn}. 
Another  approach  to  the  confidence  interval  problem  is  to  treat  the 
sequence  (f(Xn)}  as  a  stationary  process.  The  expectation  is  that 

(4.4)  n"1/2(  l  f(X.  )  -  nuf)  ->  N(0,  ol) 

k-0  1 
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where 


(4.5)  a\  -  E  (f(X  )}2  +  2  \  E^  f(XQ)  f(\)  . 

k- 1 

A  confidence  interval  procedure  can  then  be  obtained  by  using  the  CLT 

2 

(4.4),  in  conjunction  with  a  consistent  estimator  of  a^i  see  FISHMAN 
(1978)  for  a  description  of  procedures  based  on  this  idea. 

We  can  now  apply  a  theorem  from  regenerative  process  theory  to 
determine  conditions  under  which  statements  (4.4)  and  (4.5)  are 
correct.  Suppose  that  the  series  defined  by  (4.5)  converges 

-3.1.  C 

absolutely.  Then,  (4.4)  and  (4.5)  hold,  provided  Eu,  Rj  <  and 
E^|f(X)|  <  ®,  for  some  5  >  0;  see  Proposition  7.6  of  [11],  and  use 

3+  c 

a  "geometric  trial"  argument  to  prove  that  E  R.  0  <  <*>  implies  that 

*A  1 

-  3 

that  E  x_  <  ®. 

H  2 

5.  Optimal  Splitting  Measures 

Our  development  in  Sections  2  through  4  shows  that  Harris  chains 
can  be  embedded  in  a  weakly  regenerative  environment.  Furthermore,  a 
special  class  of  Harris  chains  can,  in  fact  be  made  regenerative.  We 
have  seen  that  the  regenerative  structure  can  be  used  to  advantage  in 
constructing  confidence  intervals  for  the  steady-state  simulation 
problem. 

However,  this  structure  can  be  useful  in  analyzing  a  number  of 
other  simulation  related  questions.  Procedures  that  make  heavy  use  of 
regenerative  structure  have  been  developed  for  quantile  estimation 


20 


(IGLEHART  (1976)),  extreme  values  (IGLEHART  (1977)),  sequential 
stopping  rules  for  simulation  (LAVENBERG  and  SAUER  (1977)),  and  selec¬ 
tion  of  the  "best”  system  (RUBINSTEIN  (1980)).  Given  that  a  regenera¬ 
tive  procedure  is  adopted  for  handling  a  given  problem,  a  simulator 
often  has  a  choice  between  two  or  more  sequences  of  regeneration  times 
for  the  system  under  study.  The  concensus  among  simulators  is  that 
one  should  use  the  regenerative  sequence  which  minimizes  the  expected 
time  between  regenerations  (e.g.  [10]). 

In  this  section,  we  will  therefore  consider  problems  associated 
with  construction  of  a  splitting  measure  which  gives  rise  to  the 
maximal  number  of  regenerations;  such  a  measure  will  be  called 
optimal. 

Let  {Xn}  be  a  Harris  chain  possessing  an  invariant  proba¬ 
bility  rc,  and  suppose  is  non-empty.  Then,  {(X^,  n  >  0) 

is  a  regenerative  process  with  E  x-  <  «  (Proposition  4.15,  [12]), 

H  * 

and  hence,  by  the  strong  law  for  regenerative  processes  (Theorem  3.1, 

111]) 


(5.1) 


X-1' 


/n  -*•  c 


P  a.  s. 


where  c  is  a  constant  yet  to  be  determined.  Note  that  the  quantity 
on  the  left-hand  side  of  (5.1)  is  the  number  of  regenerations  over  the 
first  n  time  units.  Applying  bounded  convergence  to  (5.1)  shows 
that  c  is  given  by 
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(5.2) 


c 


“  l  P,A  -  l)/n  . 

n  ®  k»l  ^ 


By  the  Markov  property, 

V6k  ■ 11  ■  v\_i<8  ■  »>  -  yx<Vi»  ■  y<«Vi»  • 

Hence,  by  applying  bounded  convergence  to  (3.1),  (5.2)  implies  that 

n-1 

c  -  lim  V  E  \(X._,)/n 

n  ®  k“0  ^ 

»  /  \(x)  x(dx) 

Before  proceeding,  we  point  out  that  if  (\,i)>)  satisfies  (2.2)(i) 
with  m  *  1,  and  if  /  \(x)  n(dx)  >  0,  then  there  automatically  exists 
A  such  that  (2.2) (ii)  through  (iv)  are  satisfied  (see  Proposition 
4.3,  [12]),  and  thus  (4>,\)  c  Mi . 

(5.3)  THEOREM.  There  exists  a  pair  ($,X)  e  Mj  which  maximizes 

/  X(x)  x(dx) 

over  all  pairs  ($,X)  t  Mi. 

Proof.  Since  X(x)  £  1,  it  must  be  that 

s  *  sup{J  \(x)  u(dx):  there  exists  d>  with  (4>,X)  c  M^} 
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is  less  Chan  or  equal  Co  1.  Since  Che  problem  is  Crivial  if 


is  empCy,  we  may  assume  ChaC  Chere  exisCs  a  sequence 
such  ChaC 


<  W  ‘  2, 


s(dx) 


is  increasing  Co  s.  Now,  observe  ChaC 

it(B)  -  /  P(x,B)  x(dx)  2  /  ^(x)  ^(B)  x(dx)  =  \  4^(8) 

There  exiscs  p  such  chac  for  n  2  P.  Xn  2.  s/2,  and  chus ,  for 
n  >  p, 


(5.4)  4,  (B)  <  it(B)/\  <  2it(B) / s  . 

n  —  n  — 

We  now  recall  chac  u  is  a  probabiliCy  on  Che  Borel  secs  of  a 
compleCe,  separable  meCric  space,  and  hence  it  is  a  cighc  measure 
(Theorem  1.4,  [2]).  Thus,  for  each  e  >  0,  Chere  exisCs  an  ^-compacC 
sec  such  ChaC  x(K^)  >  1-e.  LeCCing  K^_  denoCe  Che  sec  comple- 

menC  of  Kc ,  we  see,  by  (5.4)  chac 

4  (Kc)  <  2x(KC)/s  <  2e/s 
n  e  —  e 

for  n  2  P*  and  hence  <t>nC Ke )  >  1-2e/s  for  n  2  p,  proving  ChaC 
{ n »  n  2  *1  is  a  CighC  family  of  probabiliCy  measures.  Then,  by 
Prohorov's  Cheorem  ([2],  p.  35-41),  Che  family  {4>n;  n2  1)  is 
sequencially  compacC.  Hence,  Chere  exiscs  a  subsequence  {n^}  such 
ChaC 
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(5.5) 


=  >  « 


where 


is  a  probability  on  (E,£).  But 


( 4>  ,  X  )  c  M, ,  so  that 

\  “k  _1 


P(x,B)  >  (x)  6  (B) 

~  \  \ 

for  all  k  and  this  implies  that 


(5.6)  P(x,B)  >  lim  X  (x)  4>  (B)  >  lim  X  (x)  lim  <t>  (B)  , 

\  ^  “  \  -  Dk 


Let  X(x)  be 

(5.6)  implies 

(5.7) 


the  E-measurable  function  lim  X  (x) ,  and  observe  that 

\ 

that 


P(x,B)  >  X(x)  Ha 


°k 


(B) 


for  B  e  E.  In  particular,  (5.7)  holds  for  all  open  sets  G  c  _E,  and 
thus 


(5.8)  P(x,G)  >  X(x)  lim  4>  (G)  >  X(x)  41(G)  , 

-  °k  " 

the  last  inequality  by  the  weak  convergence  relation  (5.5)  (see  [2], 
p.  12).  For  a  fixed  x,  let  v(«)  be  the  signed  mesure  P(x,*)  - 
X(x)4>(*).  Applying  Hahn's  decomposition  to  v  (ROYDEN  (1968))  allows 
one  to  write 
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where  v+,  v-  are  mutually  singular  non-negative  finite  measures. 
The  measures  v+  and  v-  are  outer  regular  (see  p.  402,  3REIMAN 

(1968))  and  therefore,  given  B  c  Z,  there  exist  sequences  of  open 

12  12 

sets  G,  ,  G.  such  that  B  =  G.  ,  3  =  G. ,  and 

j  J  J  1 


v+(B)  =  lim  v+(Gj) 
v  (B)  “  lim  v  (Gj)  • 


Let 


1  2 

G.  n  G. ,  and  observe  that  B  =  G. 

j  J  J 


v+(B)  _<  lim  v+(G  )  <  lim  v+(G*)  =  v+(B) 
v'(B)  <  lim  v"(G  )  <  lim  v"(G^)  =  v~(B)  . 

Hence,  since  the  Gj  are  open,  (5.8)  implies  that 

v(B)  =  lim  v(Gj)  >_  0 

and  thus  v  is  a  non-negative  measure.  By  the  remark  preceding 
Theorem  5.3,  this  proves  that  (\,$)  c  .  On  the  other  hand,  using 

the  fact  that  \(x)  _<  1  allows  us  to  apply  Fatou's  lemma  to  prove 
that 


s  *  lim  J  \^(x)  x(dx)  <(  J  lim  X  (x)  u(dx)  »  j  X(x)  n(dx) 


concluding  the  proof  of  the  theorem.  I 
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Thus,  there  exists  a  splitting  measure  $  which  is  optimal  over 
M^.  Unfortunately,  the  above  existence  proof  is  non-cons tructive; 
furthermore,  the  measure  x  is ,  in  general,  unknown  to  a  simulator, 
and  hence  $  can  not  be  calculated  prior  to  initiating  the 
simulation.  However,  under  reasonably  general  conditions,  one  can 
find  a  sub-optimal  splitting  measure  $  which  is  optimal  in  a  certain 
"maxi— min"  sense.  Let 

?(A)  **  {4:  4  is  a  probability  on  (E,jl),  4(A)  =■  1}  for  A  c  Z  , 
and  consider  the  maxi-min  problem, 

J 

(5.9)  max  min  /  X(x)  p(dx) 

(\,6)  c  M  4  c  Y(A) 

Given  that  m  is  unknown,  this  problem  amounts  to  maximizing  with 
respect  to  the  "worst"  possible  distribution  of  x  over  A.  By 
letting  4  range  over  the  "point  mass"  probabilities  on  A,  we  see 
that  (5.9)  is  equivalent  to 

(5.10)  max  inf  \(x)  . 

(X  ,<J>)  c  x  c  A 

By  (2.2)(ii),  it  is  clear  that  there  exist  feasible  solutions  (\,i>) 
to  (5.10),  with  positive  objective  value,  for  certain  subsets  A. 

So  simplify  the  analysis  of  (5.10),  we  re-write  it  in  terms  of 
t(*)  “  X4>(*)  as  follows: 
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(5.11) 


max  t(E) 

x 

subject  to  P(x,«)  x(*)  for  x  c  A.  To  retrieve  (4>,X),  set 

\  -  t(E)  ,  $>(  •)  =  x(  •  )/\. 

(5.12)  PROPOSITION.  If  x(A)  >  0,  there  exists  a  measure  x  which 
solves  (5.11). 

Proof .  Note  that  for  any  feasible  x, 

x(‘)  n(A)  <  f  P(x,«)  u(dx)  _<  x(*) 

A 

and  hence  x(*)  is  absolutely  continuous  with  respect  to  x(0,  with 

derivative  r(*)  (say).  Now,  .E  is  the  Borel  field  of  a  separable 

metric  space,  so  tht  _E  is  generated  by  some  countable  family  of 

subsets  Bq,  B  ,  Bfl,  ...  ( DELLA CHERIE  and  MEYER  (1978),  p.  10-1 i). 

Letting  B  represent  the  o-field  generated  by  B. ,  B,,  ...,  B  , 

— n  U  1  n 

every  y  <r  E  belongs  to  a  unique  atom  A^(y)  c  J?n,  so  that 

P(x,B)  -  /  p  (x,y)  n(dy)  +  P  (x,B)  ,  for  B  e  B  , 

B  n  s  -n 

where  Pn(x,y)  =  P(x,  A^(y) )/tc (A^(y) )  (0/0  interpreted  as  zero)  and 

P  (x,*)  is  the  x-singular  component  of  P(x,*).  Let  p  (•)  be  the 
simple  function  defined  by 

p  (•  )  “  inf  {p  (x,  •)’•  x  c  A) 
n  n 


and  observe  that 
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p  (O  -  inf  P(x,  A  (.))/n(A  (•))  >  x(A  (*))/x(A  (•))  -  r  (•) 
n  .  n  n  —  n  n  n 

x  c  A 

where  rn(*)  is  Che  it-derivative  of  x  o a  Bn.  Let  p(«) 

»  lim  Pn(*)»  p(x,*)  “  lim  pn(x,*),  r(*)  -  lim  r^(>).  By  a  classical 
differentiation  theorem  (see  DOOB  (1953),  p.  612),  r( •)  -  r(*)  -z 
a.e.,  and  p(x,*)  is  the  it-derivative  of  P(x,*).  Letting 

x*(B)  -  /  p(y)  u(dy)  , 

B 

it  is  clear  that  x*(E)  x(E)  for  any  feasible  x ,  since 
p(*)  r(*).  On  the  other  hand,  for  any  x  e  A, 

/  p(y)  it(dy)  -  /  lim  inf  p  (x,y)  it(dy)  <  /  lim  p  (x,y)  it(dy) 

B  BxcA  B 

<  P(x,B) 

and  thus  x*  is  feasible.  The  measure  x*  therefore  solves  (5.  II).  II 

The  proof  of  Proposition  5.12  constructs  the  maxi-min  optimal 
measure  x*,  given  the  invariant  probability  it.  However,  the  con¬ 
struction  of  x*  can  often  be  given  in  the  absence  of  explicit  know¬ 
ledge  of  it.  Suppose,  for  example,  that  the  measures  P(x,»)  are  all 
equivalent  to  some  common  measure  q,  as  x  ranges  over  A  (two 
measures  are  equivalent  if  they  share  the  same  sets  of  measures  0). 
Then,  a  trivial  adaptation  of  the  proof  of  Proposition  5.12  shows  that 
x*  is  given  by 
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(5.13) 


r*(B)  ■  /  11m  inf  p  (x,y)  n(dy) 

B  x  c  A  n 

where  p^Cx,*)  Is  the  q-derivative  of  P(x,»)  on  B^. 

(5.14)  EXAMPLE.  Let  {Un;  n  _>  0}  be  a  sequence  of  Independent, 
Identically  distributed  real-valued  random  variables  with 


a  exp(ay)  dy  ,  y  <  0 

P{U  c  dy}  - 
n 

a  exp(-y)dy  ,  y  2  0  . 

where  dy  is  Lebesgue  measure  and  a  =  a/(<r+l).  The  process  =  x, 
»  max{Wn+Un,  0}  is  a  Markov  chain  corresponding  to  the  waiting 
time  process  of  the  M/M/1  queue  (see  IGLERART  (1971)).  If  A  »  [0,b], 
then  (5.13)  shows  that  the  Lebesgue  derivative  of  z*  is  given  by 

a  exp(a(y-b)),  0  y  £  ab 

r(y)  “ 

a  exp(-y),  y  >  ab 


We  should  point  out  that  in  actual  simulation  applications,  one  would 
use  the  maxi-min  optimal  measure  r*  as  follows.  Letting  Ajj 
represents  the  atoms  of  Bjj,  put 


(5.15) 


Pn(x) 


min 

A  c  A 
i  — n 


P(x,Ai)/-c*(Ai) 


which  we  note  is  a  decreasing  sequence  of  measurable  functions  bounded 
below  by  1.  If  we  denote  the  r*-derivative  of  P(x,»)  on  by 

Pn( x , * ) ,  then 
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P(x,B)  /  11m  p  (x,y)  x*(dy)  _>  P  (x)  x*(B)  _>  t*(B)  , 

B  n  n  - 

so  that  ($,\)  c  Mi  where 

X(x)  »  lim  p  (x)  x*(E),  4 >(•)  -  x*(‘)/x*(E) 
n 

The  algorithm  of  Section  A  could  then  be  exploited  to  obtain  a 
regenerative  sequence.  Of  course,  in  many  applications,  \( •  )  can  be 

A 

calculated  directly  from  the  densities  p(x,y)  (the  x*-derivative  of 
P(x,« ))  without  passing  through  (5.15),  namely  via 

(5.16)  X.(x)  »  inf  p(x,y)  x*(E) 

y 

The  construction  of  \(x)  via  (5.15)  is  necessary  to  deal  with  the 
case  where  one  is  working  with  a  badly  behaved  version  of  p  (with  a 
bad  version,  (5.16)  might  lead  to  \(x)  =  0). 


6.  Passage  Times  for  Harris  Chains 

In  Sections  2  through  5,  we  have  developed  simulation  methodology 
for  Harris  chains,  and  investigated  some  related  questions.  The 
theory  considered  so  far  has  pertained  exclusively  to  study  of 
problems  associated  with  simulation  output  analysis  of  functions  of 
the  form  f(Xn) ,  where  {Xn}  is  the  simulated  chain  and  f  is  a 
real-valued  ^  measurable  function.  Recently,  however,  the  "passage 
time"  problem  has  attracted  considerable  attention  in  the  literature; 
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see,  for  example,  IGLEHART  and  SHEDLER  (1980).  As  we  shall  see,  the 
quantities  of  Interest  in  the  passage  time  setting  can  not  be 
expressed  in  the  form  f(Xn) ,  and  hence  we  must  develop  some  addi¬ 
tional  theory  to  handle  this  problem. 

To  simplify  the  exposition,  we  shall  consider  only  a  special  case 
of  the  general  passage  time  problem.  However,  the  structure  to  be 
imposed  here  provides  enough  flexibility  so  as  to  cover  the  large 
majority  of  practical  cases  of  interest.  We  start  by  assuming  that 
Un>  is  a  Harris  chain  with  c-finite  invariant  measure  x.  Given 
four  E-measurable  sets  A^ ,  A2 ,  ,  B2 ,  assume  that: 


(6.1)  1)  px{\  <  <*>}  *  1  for  all  x,  where  Sq  *  0,  and 

Sk+1  -  tnfln  >  Sk:  (X^,  X„>  r  A,  *  Aj) 

il)  P  {T.  <  «}  ■  1  for  all  x,  where  T„  ■  S,  ,  and 
x  k  0  1 

Vi  ‘  l"£'»  >  V  <Vr  V  *  Bi  *  V 

111)  P  (S,  <  T,  <  S,  <  T_  <  S,  <  T,  <  •••)  -  1  for  all  x. 
x  1  1—2  2—3  3  — 

The  passage  time  problem  is  concerned  with  obtaining  estimators 

associated  with  the  amount  of  time  required  for  the  chain  (X  , ,  X  ) 

n— 1  n 

to  make  a  "passage"  from  initial  configuration  A^  x  A^  to  terminal 
configuration  x  Condition  (6.1)(iii)  guarantees  that  the 

start  times  {Sj}  and  terminal  times  (Tj)  alternate,  so  that 
passage  time  quantities  Rfj(f)  may  be  defined  via  the  formula 
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V1 

R  (f)  -  l  f(X  )  , 

J-sa  J 

for  f  a  real-valued  _E-measurable  function. 

(6.2)  EXAMPLE.  Let  {9t:  t  0)  be  a  semi-Markov  process  taking 
values  in  state  space  I  “  {0,  1,  n).  Let  P  be  a  stochastic 

matrix  which  describes  the  evolution  of  ®Tn+>  where  Tn  are  the 
successive  state  transition  times  of  the  process  {0t}.  Also,  let 
F^(«)  be  the  distribution  of  the  holding  time  in  state  i,  or 
alternatively,  let 


Fi<* 


)  -  P{  T 


n+i 


-  T  < 
n  — 


T  + 
n 


i> 


The  process  (9t)  corresponds  to  a  Markov  chain  (Zn)  on 
I  x  [0,<=°),  with  the  transition  kernel 


P((i,x),  {j}  x  [0,y}>  »  p±j  Fj(y); 

the  process  Zn  “  (sn,cn)  records  the  sequence  of  states 

visited,  together  with  the  sequence  of  consecutive  holding  times. 

Suppose  that  P  is  such  that  there  exist  four  subsets  , 

B.,  B_  on  I  such  that  (6.1)  is  satisfied,  with  the  chain  (s  } 

1  i  n 
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playing  the  role  of  {X^} ,  and  A^ ,  A2,  B ^ ,  B2  replacing  A^  A2,  B ^ , 

B».  Tlien,  the  chain  Z  ■  (s  ,c  )  itself  satisfies  (6.1),  using  the 
2  n  n  n 

sets  A^  "  A^  x  [0,®),  B^  -  B^  x  (0,®). 

Assuming  that  one  requires  the  time  required  for  { 0t }  to  make 
a  passage  from  A2  to  B2 ,  gives  that  6t  visited  Aj 
immediately  before  A2,  and  that  a  visit  to  B^  immediately 
precedes  a  visit  to  B2  (see  [18]  for  the  significance  of  A^ , 

Bp,  the  quantity  of  interest  is  R^f),  where  f(s,c)  =*  c. 

The  interest  in  passage  times  arises  from  the  fact  that  passage 
time  quantities  are  often  important  design  criteria  in  development  of 
networks  of  queues. 

We  proceed  by  first  analyzing  the  chain  Y  *  (X  ,  X  , , ) . 

n  n  rvr  1 

(6.3)  PROPOSITION.  The  chain  {Yn:  n>0}  Is  Harris  recurrent, 
provided  {Xq}  Is  Harris.  If  {Xq}  possesses  an  invariant 
probability,  then  {Yn}  does  also. 

A 

Proof.  Let  A,  X,  and  $  be  as  in  Definition  2.2.  Put  A  »  E  x  A 

and  note  that  (2 . 2) ( ill)  ensures  that  P  {Y  c  A  infinitely  often)  -  1 

pn 

2 

for  any  p  on  the  product  space  E  ,  where  P  is  the  kernel  associ¬ 
ated  with  { Yn }  given  by 

P((x1,x2),  Dl  x  D2)  ■*  ^  (*2)  P(x2,  D2)  . 

Of  course,  for  (x^,  x2)  c  A, 
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Pnri'1((x1  ,x2) ,  D1  X  d2> 

“  Px,{(V  c  D1  *  V  -  K  /  ♦(dxi)  f  P(xj  >dx2) 

*  D1  °2 

-  \iJs(D1  x  D2) 

and  hence  Yn  satisfies  (2.2)(i)  through  (iii) ,  substituting  A,  $, 
nri-1  in  place  of  A,  $ ,  m.  This  suffices  to  prove  Harris  recurrence 
(see  [12],  Section  4). 

Also,  given  that  it  is  the  unique  o“finite  invariant  measure  of 

(X  },  it  is  simple  to  verify  that  the  measure  it  given  by 
n 

(6.4)  n(D  x  D  )  =*  /  it(dx.  )  j  P(x.  ,  dx„) 

°2  °2 

A 

is  both  a-finite  and  invariant  for  P.  A  trivial  consequence  of  the 

A 

representation  (6.3)  is  that  it  is  a  probability  provided  the  same  is 
true  of  x.  B 

Let  {Xn}  be  a  Harris  chain  and  let  D  c  _E  satisfy  it(D)  >  0, 

where  it  is  the  invariant  measure  of  {X  }.  Then,  setting  T.  *  -1, 

n  u 

Vi  ■  ln£("  >  V  x„ «  D>  • 
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we  have  that  P^CE^  <  ®)  -  1  for  all  x  (see  (2.3)).  Let  {V^}  be 
the  "interblock"  chain 


S>„-  . *r  > 

n  n+1 


taking  values  in  EJ  x  D. 


(6.5)  PROPOSITION.  The  interblock  chain  {pn}  is  Harris  recur¬ 
rent,  provided  {Zn}  is  Harris.  If  {Zq}  possesses  an  Invariant 
probability,  so  does  {pn}. 


Proof.  Let  u  be  the  invariant  measure  of  {X^}.  By  (2.3),  ti(C)  >  0, 
implies  that 


(6.6)  Px^Xn  c  C  finitely  often)  »  1  for  all  x  . 


For  F  a  measurable  subset  of 


x  D, 


let 


**(F)  -  c  D,  (Xj,  ...,  XJc)  «  F,  Tj  -  k)  . 

Assuming  that  7t*(F)  ^  0.  choose  6  sufficiently  small  that 


CF  -  {x  e  D:  Px{ (Xj . Xfc)  c  F;  Tj  -  k)  >  5} 
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has  positive  ir-measure.  Of  course,  by  (6.6)  Xj>Q  visits  Cp 
infinitely  often  Px  a.s.  Let  Uq  “  0, 


Vl  '  lnf<J  >  V*!  XJ  6  CF> 


and  observe  that  the  Uj's  are  a  subsequence  of  the  T^'s,  so  that 


Px{pk  t  F  for  all  k  >  1} 


<  P{(X„  WF  for  l<k<n) 

x  uk  1  k+1 


where  V^+1  =  inf{n  >  0^:  <r  D) .  Then,  by  the  strong  Markov 

property  applied  at  time  Un,  the  above  term  can  be  re-expressed  as 


{(X, 


XT  ) 
A1 


l  F} ;  (Xy 


k+1 


*V  > 
k+1 


k  F 


1  <  k  <  n) 


<  (1-6)  P((X_  X_  )  i  F,  1  <  k  <  n}  , 

x  uk  1  k+1 

the  1-6  term  resulting  from  the  fact  that  Xun  c  Cp.  Repeating  this 
argument  (n-1)  more  times  shows  that 


Px{Pk  i  F  for  all  k  >  1}  <  (l-6)n 

and  hence  visits  F  Px  a.s.  for  all  x.  This  implies  that 

(Pq}  satisfies  (2.3)  with  it*  playing  the  role  of  v,  proving  Harris 
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recurrence.  Finally,  given  that  n  is  an  invariant  probability  for 
{Xq},  it  is  trivial  to  verify  that 

(6.7)  **(•)  -  P^{X0  c  D,  iXl . XT  )  e  •} 

is  an  invariant  probability  for  {pn}.  5 

Returning  now  to  the  passage  time  problem,  we  can  combine 
Propositions  6.3  and  6.5  to  obtain  the  following  result. 


(6.8)  THEOREM.  Suppose  that  {Xq}  is  a  Harris  chain.  Then,  if 
{ Sfc)  is  defined  through  (6.1),  the  chain 


) 


is  a  Harris  chain.  Furthermore,  if  (Xq)  possesses  an  invariant 
probability,  then  so  does  (5Q}. 

Because  of  the  alternation  condition  (6.1)(iii),  it  follows  that 


R”(£>  ’  J  £<x'‘> 

n 

can  be  taken  to  be  a  measurable  function  of  the  {Yn}  process,  for 
any  E-measurable  f.  Since  is  a  functional  of  a  Harris 

chain,  one  an  immediately  apply  the  general  theory  of  Harris  chains  to 
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obtain  sufficient  conditions  under  which  strong  laws  and  central  limit 
theorems  hold  for  partial  sum  processes  associated  with  {g(Rn(f))}, 
for  real-valued  Borel  measurable  g  (typical  choices  of  g  might  be 
g(x)  *  x^).  To  be  precise,  it  is  evident  that  if  x  is  a 
probability,  then 

n 

(6.9)  £  g(Rk(f))/n  -  Eit(g(R1(f))} 

Px  a.s.  for  all  x  c  E,  provided  that  the  right-hand  side  of  (6.9) 
is  finite;  the  identification  of  the  limit  in  (6.9)  follows  from 
(6.7).  If  we  assume  a  further  moment  condition,  then  the  weak  regen¬ 
erative  embedding  of  Section  2,  together  with  the  CLT  of  Section  3, 
proves  that  there  exists  a  constant  such  that 


n1/2[  l  g(5L(f))  -  E  g(R  (f))]  *»>  N(0,  a2)  . 
k-1  *  xi 

We  can  also  obtain  a  simple  sufficient  condition  for  finiteness 
of  g(Ri(f)),  in  the  case  where  g(x)  =»  x. 

(6.10)  PROPOSITION,  If  v | f |  <  •,  then  ) |  <  •. 

Proof.  Let  h(X^_^,  X^)  -  |f(X^_^)|  +  |f(X^)|  and  observe  that  if  x 

is  the  invariant  measure  of  Y  ■  (X  ,,  X  )  then  by  (6.4), 

n  n- 1  i 
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/  h(x1,x2)  uCdXp  dx2)  -  E^  h(X^_1 ,  X^)  -  2-a  j  f  j  <  ® 

so  that  h  is  n-integrable.  Now,  the  sequence  of  start  times  { S^} 
constitutes  the  consecutive  hitting  times  of  (A^,A2)  for  the  Harris 
chain  Yn.  A  well-known  construction  of  the  invariant  measure  of  a 
Harris  chain  (see  [20]  ,  p.  32)  shows  that 

Sk+1 

*(D)  =  E.{  V  I  (Y  )}/u(A.  x  A  )  . 

*K+l 

A  classical  approximation  argument  then  proves  that 

S*+l 

(6.11)  /  h  it(dy)  -  Ea{  l  h(X  X  ) }/n(A  *  A  ) 

r.  j=Sk+l  J  J 

and  hence  if  n | f j  <  “,  the  right-hand  side  of  (6.11)  is  finite.  But 

Sk+1  Sk+1 

E  Rn(f)  <  E-{  I  |f(X  )|}<E.{  )  h(X  X  )}  <  «  , 

j=Sk+1  3  it  j“Sk+1  J  J 

proving  our  assertion.  I! 


7.  A  Nonlinear  Storage  Model 

In  this  section,  we  study  a  discrete  time  storage  model  which  is  a 
vector  version  of  a  continuous  time  process  studied  by  (JINLAR  and 
PINSKY  (1971),  and  by  HARRISON  and  RESNICK  (1976).  We  shall  prove 
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that  the  model  gives  rise  to  a  Harris  chain  possessing  an  invariant 

probability;  the  argument  will  Illustrate  some  techniques  applicable 

to  proving  Harris  recurrence  in  applied  probability  models. 

The  model  that  we  shall  consider  concerns  a  finite  family  of 

interconnected  reservoirs.  Suppose  that  S  “(SO),  ....  S  (!)) 

n  n  n 

are  the  storage  levels  at  time  n  in  each  of  the  X  reservoirs. 

Over  the  interval  (n,  n+1],  the  l  reservoirs  receive  additional 

content,  from  an  external  source,  in  the  amount  ZjH-]..  Given  the 

storage  levels  Sn  at  time  n,  a  decision  is  made  to  release 

amounts  R  -  (R  (1),  ....  R  (X ) )  over  the  time  interval  (n,  n+1]. 
n  n  n 

Assuming  that  the  release  rule  R  “  F(S  ),  the  above  formulation 

n  n 

leads  to  a  recursion 


(7.1) 


n+1 


S  +  Z 


n+1 


-  r(sa) 


We  now  make  the  following  assumptions  on  F: 


(7.2)  i)  F  is  strictly  increasing, 

ii)  s-r(s)  is  strictly  increasing, 

iii)  F  is  continuously  differentiable, 

iv)  T(0)  -  0. 

By  strictly  increasing,  we  mean  that  is  positive,  where  F^ 

is  the  ith  component  of  F» 

We  introduce  a  stochastic  element  into  the  model  by  assuming  that 
(Z^:  n  1}  is  a  sequence  of  non-negative  independent  and  Identically 
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distributed  random  vectors.  The  following  requirements  are  imposed  on 
the  joint  distribution  function  F(z)  of  { Zn } t 

(7.3)  i)  E2  <  lim  F(s), 

•  n  - 

3  s5  -*•  ® 

ii)  F  has  a  continuous  Lebesgue  density  component  f  which 
is  positive  on  Be  =  {x:  F(a-e)  <  T(x)  <  T(a+e)}  for 
some  z  >  0,  where  F(a)  >  EZn  (inequality  here  is 
componentwise) . 

As  previously  mentioned,  this  model  bears  close  resemblance  to 
some  recently  studied  continuous-time  storage  processes;  it  is  also 
closely  related  to  the  one-variable  discrete-time  storage  model  of 
BATHER  (1962). 

(7.4)  THEOREM.  Under  conditions  (7.2)  and  (7.3),  {Sq}  is  a  Harris 
chain  with  invariant  probability. 

Proof.  For  q  >  0,  let  A  -  (s:  F(s)  <  EZ  +n}  “  (s:  s  <  F  *(EZ  -rn)}  , 
-  T)  —  n  —  n 

which,  by  (7.2)(iii)  and  (7.3)(i),  is  non-empty  and  compact  for  r| 
sufficiently  small.  Now,  let  g(s)  »  Dsn  (3»3  -  sum  of  absolute 
values  of  components)  and  observe  that 
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F 


1 


Eg  g(SL)  -  EJs  +  Zx  -  r(s)J 


l 

-  l  E{s(i)  +  Z  (i)  -  r.(s)} 
i-1  1  l 


l  l 

-  I  S(i)  +  l  (ez  (i)  -  r,(s))  <  BSB  -  Jin  -  g(s)  -  Jin 
i-1  i-1  1  1 


for  s  l  A^,  where  the  second  equality  follows  by  non-negativity  of 
s-r(s).  Also, 


sup  ec  g(s.)  -  sup  ns-r(s)B  +  n ez  n  <  ® 

s  c  A  1  s  c  A  n 

h  h 

by  compactness  of  A^  ,  so  that  by  Theorem  6.1,  TWEEDIE  (1976),  the 
set  A^  is  uniformly  positive.  Hence, 


(7.5)  sup  E  T  <  ® 

A  8  A 

n 

n 

where  -  inf{n  >0:  S„  c  A^} .  By  (7.3)(ii),  we  can  choose  n 

and  y  <  t  sufficiently  small  that  A  n  B  »  we  will  show  that 

n  y 

there  exists  n,  5  >  0  such  that 


(7.6) 


sup 
s  c  A 


P  (T_ 
sl  B 


>  n}  <  1-6 


Let  4>(s)  -  s  -  T(s)  +  r(a~Y/2),  and  observe  that  <S  is  a  strictly 
monotone  transformation,  so  that  <£n(s)  increases  to  a  -  y/2  for 
s  c  A^.  In  fact,  <*“(0)  <  <t>n( s )  so  that  it  is  possible  to 
choose  an  n  for  which 
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3>n(s)  >  a  -  y/2 


uniformly  over  s  c  A^.  Hence,  for  s  c  A,,, 

VTB  <n>  >  Ps<Snc  V  Zl*  •* 

Y 

-  pstr(Sn)  >  r(a  -  y/2)  ;  Zj,  . 


Zn  c  By/2} 


C  BY/2} 


-  Ps{$n(s)  >  a  -  y/2;  Zn  c  B  /2) 


’  W  —•  Zn  c  BY/2>  “  6  >  0  * 

proving  (7.6).  Relations  (7.5)  and  (7.6)  together  show  that  for  any 

Y  >  0, 


(7.7)  Ps{Tb  <  •}  «  1 

Y 

for  all  s.  Now,  use  (7.3)(ii)  and  note  that  for  y  sufficiently 
small 


(7.8)  Ps(S1  c  C}  -  P{Z1  +  s  -  r(s)  e  C}  >  X  /  dy 

c  "\n 

for  all  s  e  By,  where  X  is  the  lower  bound  on  the  Lebesgue 
density  of  F  over  BE .  Together,  (7.7)  and  (7.8)  suffice  for 
Harris  recurrence,  in  light  of  (2.2). 

The  fact  that  the  invariant  measure  is  a  probability  will  follow 
as  a  consequence  of  being  compact  and  uniformly  positive. 

First,  (Sn)  is  a  Harris  chain  with  Feller  kernel;  i.e. ,  sn  s 
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implies  P(sn,«)  *—>  P(s,*)  (this  is  a  consequence  of  the  fact  that 
(7.1)  defines  Sn+j  as  a  continuous  function  of  Sn).  We  may 
therefore  apply  Lemma  5.1  of  [24],  which  shows  that  A^  qualifies 
as  a  positive  status  set.  Application  of  Proposition  4.3  of  [24] 
completes  the  proof  of  the  theorem.  I 

The  theory  of  Section  4  proves  that  under  conditions  (7.2)  and 
(7.3),  {S^}  may  be  embedded  in  a  regenerative  environment.  As  a 
special  case,  consider  the  one  variable  model  in  which  T(s)  »  as  for 
0  <  a  <  1.  Then,  {Sjj}  is  an  autoregressive  process  of  order  1 
and 

”-9>  Sn+  Vl  • 

Given  that  F(dz)  «  f(z)dz,  where  f  is  a  positive  continuous 
density,  such  an  autoregressive  process  may  be  regarded  as 
regenerative  when  appropriately  embedded.  A  natural  question  to  ask 
is  whether  {Sn}  can  be  made  regenerative  without  the  embedding. 

More  precisely,  is  the  process  { Sn}  strongly  regenerative  in  the 
sense  that  there  exist  regeneration  times  that  are  measurable 
functions  of  the  process  {SQ}  alone?  By  appealing  to  Theorem  8.17 
of  [11],  we  see  that  {Sn}  can  be  strongly  regenerative  only  if 
there  exist  sets  ,  C2  of  positive  Lebesgue  measure  such  that 

(7.10)  Pg{S1  c  dy}  *>  f(y  -  (l-a)s)dy  -  cg  g(y)  dy 


44 


for  some  function  cs,  g,  where  (7.10)  must  hold  for  all  s  £  Cj, 
y  c  C2.  Since  the  factorization  (7.10)  is  not  generally  valid,  we 
have  proven  existence  of  a  chain  that  can  be  made  regenerative  only 
via  an  embedding. 

It  is  also  worth  pointing  out  that  some  sort  of  density  assump¬ 
tion  for  {Zn}  is  necessary,  in  order  to  obtain  Harris  recurrence. 

In  Section  6  of  [12]  ,  it  is  shown  that  the  autoregressive  process 
defined  by  (7.9)  can  not  be  embedded  in  a  weakly  regenerative  environ¬ 
ment,  even  in  the  presence  of  the  condition  EZn  <  ®,  if  {Zn}  has 
a  certain  atomic  distribution. 
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